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Path- Integral- Monte- Carlo simulation has been used to calculate the proper- 
ties of a two-dimensional (2D) interacting Bose system. The bosons interact 
with hard-core potentials and are confined to a harmonic trap. Results for 
the density profiles, the condensate fraction, and the superfluid density are 
presented. By comparing with the ideal gas we easily observe the effects of 
finite size and the depletion of the condensate because of interactions. The 
system is known to have no phase transition to a Bose-Einstein condensation 
in 2D, but the finite system shows that a significant fraction of the particles 
are in the lowest state at low temperatures. 
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1. INTRODUCTION 

Recent experiments on alkali atoms cooled by laser methods and evap- 
oration in a magnetic trap have allowed the observation of Bose-Einstein 
condensation (BEC) in a harmonic potential formed by an external mag- 
netic field. By varying the trapping field so that it is very narrow in one 
dimension, it is possible to separate the single-particle states in the oscilla- 
tor potential into well-defined bands. By occupying only states in the lowest 
band one has an effective two-dimensional system. In contrast to such a 
quasi-2D system, the system considered here is genuinely two dimensional. 

In this paper we will investigate the behaviour of harmonic Bose sys- 
tems in 2D using the very powerful finite-temperature Path-Integral Monte 
Carlo (PIMC) simulation technique.^ The PIMC technique is in principle 
capable of describing systems of arbitrary interaction strength and density 
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and allows one to study the static properties of the condensed gases. The 
only fundamental uncertainty arises from the choice of the interaction poten- 
tial. Here a hard-core potential appropriate for the s-wave scattering length 
of 87 Rb was chosen. In this case the hard-core parameter is ao = 0.0043 in 
the dimensionless units of Ref. 2. The analogous three dimensional case has 
been studied previously.E 

The Hamiltonian for the system under consideration is 

N 2 N 

H = E £ + E v ^ - + y Efe* + (i) 

i=l i,j i=l 

The density matrix for N Bose particles at inverse temperature (3 = 
1 /ki,T can be written as a convolution with M intermediate density matrices 
or "time slices" at inverse temperature r = (3/M. Then the probability 
density for finding a many-particle configuration R = (n, • • • , rjy) is 

p(R,R;/3) = ^yE /••• / />(RRi^) x 

x p( Ri, R 2 ; r) • • • /)(R M -i,R p ; r)dRi<iR2 • • • dR M -i 



where R^ denotes a vector with permuted particle labels. The Metropolis 
algorithm is used to sample from this distribution. With larger M or corre- 
spondingly larger temperatures for each time slice, the intermediate density 
matrices approach the classical limit which leads to the primitive approxi- 
mation. In order to make the computation for many particles feasible, M 
can be reduced by several orders of magnitude by calculating the density 
matrix p2 for the interaction involving just two particles and approximating 
each time slice bycl 

/-rj -rW \ T~T / , \ TT fa » r i j r 'i » r j j T ) 

p(R, R ; r) = ^ ^) E[ ^(r^^r^^y.r) " (2) 



As long as only two particles interact, this is equivalent to using the primitive 
approximation. 

To calculate p2 , we note that in the harmonic potential the Hamiltonian 
for two particles decouples into a center of mass and a relative motion term 
with the latter given by 

H rel = ^ + V{r)+^y x +u>y y ) (3) 
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where p = m/2 is the reduced mass. With the Trotter breakup, the quotient 
in (^), which plays the role of a correction term, can be transformed to 

, ,.X(r)X(r') 
pHc(r,r') ' ' 4 

PiA r > r ) 

with X(r) = exp(— r//(o;^r^ + u^r^)/4) and the relative coordinates r = 
Ti — Tj and r' = r^ — r'-. pHc( r i r> ) is the density matrix for the Hamiltonian 
Hue (Eq. (||)) describing the interaction between two interacting particles 
without confining potential and pi )t i{r,r') is the density matrix for a single 
particle of reduced mass p in the harmonic potentialJa 

Calculating pHc{ r , r ') f° r a hard-core potential using an eigenfunction 
expansion yields 

PHc(ry-,r)= Y, dkR kl (r)R kl (r')e- T ^t (5) 

7 ^ .. J 

1 = — OO 

with the radial wave functions 

Rki(r) = Vk{cos5 M Ji(kr) - sm5 k iNi(kr)) and taxi5 H = ^ 

N z (fca ) 

where J; and N; are the Bessel and Neumann functions and <p is the angle 
between r and r'. 

Since the numerical evaluation of this function is quite time consuming, 
values distributed over a mesh with parameters |r|, \ r'\ and <p are computed 
for each value of r. A simple linear interpolation method proved to be 
sufficient to reach the required degree of accuracy for evaluation in the sim- 
ulation. The full correction factor (Eq. (||)) with dependence on the confin- 
ing potential can then be calculated very efficiently during the Monte-Carlo 
simulation. 

To further speed up the calculation, a boxing algorithm! is used that 
divides the space of the system into boxes with a size of at least the "healing 
length" of the pair interaction. When we compute the effects of the inter- 
action of a given particle, it is necessary to consider only a small number 
of interactions with particles in boxes neighbouring that of the particle in 
question. Furthermore, attempting only permutations with particles from 
the same box, increases the acceptance probability for these moves and re- 
sults in a more effective sampling of permutation space.! 

The value of r has to be chosen carefully for each particle density to 
accommodate both the Trotter breakup and the approximation of pair inter- 
actions (Eq. (|2|)). For the density used, tests with values for r ranging over 
orders of magnitude have been performed, showing that a further decrease 
below t = 0.01 yield the same results within statistical errors. 
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2. DENSITY PROFILES AND CONDENSATE FRACTION 

Ji 2D, with interactions, there is no phase transition to a Bose condensed 
even in a trap, nevertheless there should be a macroscopic fraction of 
particles in the ground state at finite temperature when the particle number 
is finite. 

We tentatively assume that, in 2D, the condensation fraction can be 
determined, as in 3D,u by observing a macroscopic number of particles with 
long exchange cycles in the lowest state of their subsystem. For each tem- 
perature a characteristic length of permutation cycles Iq can be chosen so 
that the density profiles for particles on cycles longer than Iq are essentially 
the same. Particles on shorter cycles will have a broader density profile. 
Particles on cycles longer than Iq are identified with the condensate (Fig. 1). 
What is plotted is the average distribution in one coordinate after integrat- 
ing over the other coordinate. We have fit the lowest temperature curve in 
Fig. 1 with the infinite-N solution of the zero-temperature GP equation.! 
The interaction coefficient has been determined by the fit since, in 2D, there 
is no straightforward connection between a hard-core interaction and a con- 
tact pseudopotential. We observe no obvious bimodality, due to separate 
distributions of condensed and non-condensed particles in the overall den- 
sity in Fig. 1; if a kink did exist in the density it would likely be washed out 
by the integration over one coordinate. 

In two dimensions interactions seem to lead to a comparatively stronger 
depletion of the condensate (Fig. 2) than in three dimensions (Cf., Ref. 3). 
This is expected, since there is no condensate at finite temperature in the 
thermodynamic limit for the interacting system and therefore the critical 
temperature has to approach zero (for large numbers of particles) when 
interactions are turned onU 



3. SUPERFLUID FRACTION 

Recently one of the authorsi showed that the Hartree-Fock-Bogoliubov 
equations indicated that there is a phase transition in 2D, although it cannot 
be to the BEC state. Perhaps this transition is of the Kosterlitz-Thouless 
(KT) typei and involves superfluidity. The superfluid fraction can be related 
to the mean square surface area enclosed by Feynman paths. Sindzingre 
et. al. have shown thara 

Am 2 {A 2 ) 

with A the area swept out by the paths and I c the classical moment of 
inertia. The average is taken over configurations in the simulation. 
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The resulting superfluid densities are very small (Fig. 2). A smaller p s 
than in the translationally invariant system is expected because the forma- 
tion of the superfluid begins in the middle of the potential well where the 
contribution to the moment of inertia is small. Paths with different orien- 
tations will contribute with different signs to the area making cancellation 
possible. 

The simulation results for the 2D system without confining potential!^ 
are in agreement with the KT theory for the superfluid transition. If this 
description is also appropriate with confining potential, the vortex picture 
of the KT-transition suggests an additional mechanism leading to a decrease 
of pjj : Superfluidity is destroyed by dissipation through vortices which are 
not paired. At low temperatures vortices form pairs which unbind at higher 
temperatures. In the non-uniform system the two vortices forming a pair 
will in general experience a slightly different potential leading to imperfect 
pairing and thereby a lowering of the superfluid fraction. 
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Fig. 1. Density profiles for 1000 particles at temperatures T/T c = 
0.84,0.56,0.42,0.28,0.21 for u> = 0.15. The condensate parts are displayed 
on the left and the profiles for all particles on the right. All profiles are 
normalised to unity. Temperatures are given with reference to the criti- 
cal temperature T c of the ideal system in the thermodynamic limit. The 
smooth dotted curve through the 0.21 data is the infinite-N solution of the 
GP equation with interaction strength determined by the fit. 
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Fig. 2. Left figure: Condensate fraction for the interacting system with 
N=1000 and u = 0.15 and comparison with theoretical prediction for the 
ideal gas in the thermodynamic limit (TDL) and with finite size corrections 
for 1000 particles. Right figure: Dependence of the superfluid fraction on 
temperature for 1000 particles. 
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